Running univariate and bivariate LAVA
Aaron Wagen & Regina H. Reynolds
UCLAim: run univariate and bivariate tests
1 Background
Neurodegenerative diseases are often defined by the predominant neuropathology (e.g. synucleinopathies by \(\alpha\)-synuclein deposition), but more often than not, they also present with co-pathologies, suggesting that they might share common pathogenic pathways (PMID: 29878075, 30258235). This notion is supported by genetic studies, which have (i) identified shared risk loci across neurodegenerative disorders, such as APOE and BIN1 in Alzheimer’s disease (AD) and Lewy body dementia (LBD), or GBA, SNCA, TMEM175 in Parkinson’s disease (PD) and LBD and (ii) demonstrated that genetic risk scores derived from one neurodegenerative disorder can predict risk of another, as with AD and PD predicting risk of LBD (PMID: 31701892, 30953760, 33589841).
From a clinical perspective, neurodegenerative diseases are often also defined in terms of their predominant cognitive and/or motor symptoms (e.g. Alzheimer's disease by cognitive and memory impairment or Parkinson's disease by parkinsonism), but in reality, present as highly heterogenous diseases, with symptoms spanning multiple domains, including neurospychiatric symptoms (PMID: 27188934, 28332488. Indeed, a higher prevalence of depression has been observed in individuals with dementia compared to those without dementia (PMID: 30536144). Notably, a similar (albeit reversed) phenomenon has been observed in some neuropsychiatric disorders, with a higher risk of dementia diagnoses observed in individuals with schizophrenia (SCZ) versus individuals without a history of serious mental illness (PMID: 33688938, 26444987).
These clinical and/or pathological overlaps are not, however, always reflected in terms of global genetic correlations, with a recent study of global genetic correlation between neurological phenotypes demonstrating limited overlap between individual neurodegenerative disorders as well as between neurodegenerative and neuropsychiatric disorders (PMID: 29930110). One explanation for this is that global studies only consider an average across the entire genome, and thus, may fail to detect correlations confined to particular genomic regions, or where changes are balanced out over the genome. This limitation has recently been addressed with the development of Local analysis of [co]variant annotation (LAVA), which is able to evaluate local heritability over multiple traits of interest (using GWAS summary statistics) and detect local regions of shared genetic association (https://www.biorxiv.org/content/10.1101/2020.12.31.424652v2). Here, we apply LAVA to genome-wide association studies (GWASs) derived from three neurodenerative disorders (AD, PMID: 30617256;LBD, PMID: 33589841; and PD, PMID: 31701892), and three neuropsychiatric disorders (bipolar disorder, BIP, PMID: 34002096; major depressive disorder, MDD, PMID: 30718901; and SCZ, PMID: 29483656).
2 Supplementary code
Following section includes any intermediary code used in this .Rmd.
2.1 Run univariate and bivariate tests
This can either be sourced directly:
source(here::here("scripts", "02_run_univar_bivar_test.R"))Alternatively, for many phenotypes/loci it can be worth running this script directly on the cluster, using nohup, as performed below:
# Have to navigate to root project folder for script to work (as it uses here package)
cd /home/rreynolds/misc_projects/neurodegen-psych-local-corr
nohup Rscript \
/home/rreynolds/misc_projects/neurodegen-psych-local-corr/scripts/02_run_univar_bivar_test.R \
&>/home/rreynolds/misc_projects/neurodegen-psych-local-corr/logs/02_run_univar_bivar_test.log&2.2 Loading results
Once run, results can be loaded using the following code chunk:
files <-
list.files(
path =
here::here("results",
"02_univar_bivar_test"),
pattern = ".lava.rds",
full.names = T
)
results <-
setNames(
object = files %>%
lapply(., function(file){
list <-
file %>%
readRDS()
list %>%
purrr::discard(is.null) %>%
qdapTools::list_df2df() %>%
dplyr::select(-X1)
}),
nm = files %>%
basename() %>%
str_remove(".lava.rds") %>%
str_remove(".*:") %>%
str_remove(".*\\.")
)global <-
read_delim(
file = here::here("results", "01_input_prep", "ldsc_corr", "ldsc_correlations.txt"),
delim = "\t"
) 2.3 Filtering for pleiotropic genes
ref <- import("/data/references/ensembl/gtf_gff3/v87/Homo_sapiens.GRCh37.87.gtf")
ref <- ref %>% keepSeqlevels(c(1:22), pruning.mode = "coarse")
ref <- ref[ref$type == "gene"]
# genes to check
genes <-
c("APOE", "TMEM175", "SNCA", "BIN1", "GRN", "KIF5A")
# construct granges object with LD blocks for overlap with reference gtf
loci_gr <-
loci %>%
GenomicRanges::makeGRangesFromDataFrame(
.,
keep.extra.columns = TRUE,
ignore.strand = TRUE,
seqinfo = NULL,
seqnames.field = "CHR",
start.field = "START",
end.field = "STOP"
)
overlap <-
GenomicRanges::findOverlaps(loci_gr, ref) %>%
tibble::as_tibble()
loci_genes <-
tibble::tibble(
locus = loci_gr[overlap$queryHits]$LOC,
chr = loci_gr[overlap$queryHits] %>% GenomeInfoDb::seqnames() %>% as.character(),
locus_start = loci_gr[overlap$queryHits] %>% BiocGenerics::start(),
locus_end = loci_gr[overlap$queryHits] %>% BiocGenerics::end(),
gene_id = ref[overlap$subjectHits]$gene_id,
gene_name =
ref[overlap$subjectHits]$gene_name %>%
stringr::str_replace_all("-", "_") %>%
stringr::str_replace_all("\\.", "_"),
gene_biotype = ref[overlap$subjectHits]$gene_biotype %>% as.factor(),
gene_strand = ref[overlap$subjectHits] %>% BiocGenerics::strand() %>% as.character(),
gene_start = ref[overlap$subjectHits] %>% BiocGenerics::start(),
gene_end = ref[overlap$subjectHits] %>% BiocGenerics::end()
)
# filter df with ld block & gene info for genes of interest
loci_genes_filtered <-
loci_genes %>%
dplyr::filter(
gene_name %in% genes
)2.4 Loading plotting functions
source(here::here("R", "plots.R"))3 Methods
3.1 Running univariate and bivariate tests
The detection of valid and interpretable local genetic corrlation (\(r_{g}\)) requires the presence of sufficient local genetic signal. For this reason, a univariate test was performed as a filtering step for bivariate local \(r_{g}\) analyses. Bivariate local correlations were only performed for pairs of traits which both exhibited a significant univariate local genetic signal (p < 0.05/300, where the denominator represents the total number of tested loci, i.e. LD blocks). This resulted in a total of 1603 bivariate tests spanning 275 distinct LD blocks.
4 Results
4.1 Tabular results
4.1.1 Global correlations
print("Global correlation column descriptions:")## [1] "Global correlation column descriptions:"tibble(
column = colnames(global),
description =
c(
"Phenotype 1",
"Phenotype 2",
"The estimated genetic correlation",
"The bootstrap standard error of the genetic correlation estimate",
"The bootstrap standard error of the genetic correlation estimate",
"P-value for genetic correlation",
"Estimated snp-heritability of the second phenotype ",
"Standard error of h2 for phenotype 2",
"Single-trait LD score regression intercept for phenotype 2",
"Standard error for single-trait LD score regression intercept for phenotype 2",
"Estimated genetic covariance between p1 and p2",
"Bootstrap standard error of cross-trait LD score regression intercept"
)
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Global correlation results (p < 0.05/n_tests):")## [1] "Global correlation results (p < 0.05/n_tests):"combn <-
global$p1 %>%
unique() %>%
combn(m = 2) %>%
t() %>%
as_tibble() %>%
dplyr::rename(
p1 = V1,
p2 = V2
)
global_signif <-
global %>%
dplyr::filter(
!p1==p2
) %>%
dplyr::inner_join(combn) %>%
dplyr::filter(
p < 0.05/nrow(combn)
)
global_signif %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.1.2 Univariate results
print("Univariate column descriptions:")## [1] "Univariate column descriptions:"tibble(
column = colnames(results$univ),
description =
c(
"LD block ID",
"LD block chromosome",
"LD block start base pair",
"LD block end base pair",
"The number of SNPs within the LD block",
"The number of PCs retained within the LD block",
"Analysed phenotype",
"Observed local heritability",
"P-value from the univariate test (F-test for continuous, Chi-sq for binary)"
)
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Univariate results (p < 0.05/n_ld_blocks):")## [1] "Univariate results (p < 0.05/n_ld_blocks):"results$univ <-
results$univ %>%
dplyr::mutate(
phen = phen %>%
str_replace_all("[:digit:]", "") %>%
str_remove("\\..*") %>%
fct_relevel(
fct_disease
)
)
results$univ %>%
dplyr::filter(p < 0.05/nrow(loci)) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.1.3 Bivariate results
print("Bivariate column descriptions:")## [1] "Bivariate column descriptions:"tibble(
column = colnames(results$bivar),
description =
c(
"LD block ID",
"LD block chromosome",
"LD block start base pair",
"LD block end base pair",
"The number of SNPs within the LD block",
"The number of PCs retained within the LD block",
"Phenotype 1",
"Phenotype 2",
"Standardised coefficient for the local genetic correlation",
"Lower 95% confidence estimate for rho",
"Upper 95% confidence estimate for rho",
"Equivalent of taking the square of rho. Denotes the proportion of variance in genetic signal for phen1 explained by phen2 (and vice versa)",
"Lower 95% confidence estimate for r2",
"Upper 95% confidence estimate for r2",
"Simulation p-values for the local genetic correlation"
)
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Bivariate results (p < 0.05/n_bivar_tests):")## [1] "Bivariate results (p < 0.05/n_bivar_tests):"results$bivar <-
results$bivar %>%
dplyr::mutate(
phen1 = phen1 %>%
str_replace_all("[:digit:]", "") %>%
str_remove("\\..*"),
phen2 = phen2 %>%
str_replace_all("[:digit:]", "") %>%
str_remove("\\..*"),
fdr = p.adjust(p, method = "fdr")
)
bivar_thres <- 0.05/nrow(results$bivar)
bivar_bonf <-
results$bivar %>%
dplyr::filter(p < bivar_thres)
bivar_bonf %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.2 Global correlations versus bivariate tests
4.2.1 Text
- Global cross-trait correlations revealed 4 significant global \(r_{g}\)'s (Figure 4.2). These global \(r_{g}\)'s are very much in line with what we expect from previous literature e.g. significant positive correlations between the three neuropsychiatric disorders (PMID: 29930110). Between neurodegenerative disorders, we only observe a significant postive \(r_{g}\) between LBD and PD.
- With a Bonferroni cut-off of 0.000031, we detect a total of 77 significant bivariate local \(r_{g}\)'s across 59 distinct LD blocks, of which 11 were associated with more than one trait pair.
- The phenotype with the highest number of significant bivariate local \(r_{g}\)'s was SCZ followed by BIP. Notably, the number of significant bivariate local \(r_{g}\)'s did not appear to correlate with GWAS sample size (Figure 4.3).
- As seen from the table below and Figure 4.4, the phenotype pairs with the most local \(r_{g}\)'s were (i) BIP and SCZ and (ii) AD and PD.
- For some of these local \(r_{g}\)'s, the 95% confidence intervals (CI’s) for the explained variance (r2) included 1 (see table below), which is consistent with a scenario where the genetic signal of those pairs of phenotypes in these LD blocks is completely shared.
# proportion of correlations where CI 97.5% for explained variance (r2) = 1
# consistent with scenario where genetic signal of those pairs of phenotypes in these loci is completely shared
n_upper_r2_ci <-
bivar_bonf %>%
dplyr::filter(
r2.upper == 1
) %>%
dplyr::count(phen1, phen2)
phen_count <-
bivar_bonf %>%
dplyr::count(phen1, phen2, name = "total_n") %>%
dplyr::left_join(n_upper_r2_ci) %>%
dplyr::mutate(
prop_upper_r2_ci_equals_one =
n/total_n,
prop_upper_r2_ci_equals_one =
case_when(
is.na(prop_upper_r2_ci_equals_one) ~ 0,
TRUE ~ prop_upper_r2_ci_equals_one
)
) %>%
dplyr::select(contains("phen"), total_n, prop_upper_r2_ci_equals_one)
phen_count- Within neuropsychiatric disorders, all significant bivariate local \(r_{g}\)'s were positive (matching the overall positive global \(r_{g}\)), while within neurodegenerative disorders there was a mix of positive and negative local \(r_{g}\)'s (Figure 4.4).
- All three neurodegenerative disorders shared one significant local \(r_{g}\) with SCZ, and also PD shared one negative local \(r_{g}\) with BIP.
4.2.2 Figures
data_to_plot <-
read_delim(
file = file.path(here::here("results", "01_input_prep"),
"input.info.txt"),
delim = "\t"
) %>%
dplyr::select(-filename) %>%
tidyr::pivot_longer(
cols = cases:controls,
names_to = "group",
values_to = "total"
) %>%
dplyr::mutate(
phenotype_group =
case_when(
phenotype %in% c("AD2019", "LBD2020", "PD2019.meta5.ex23andMe") ~ "Neurodegenerative",
TRUE ~ "Neuropsychiatric"
),
phenotype =
case_when(
phenotype == "AD2019" ~ "Alzheimer's disease\n(AD)",
phenotype == "LBD2020" ~ "Lewy body dementia\n(LBD)",
phenotype == "PD2019.meta5.ex23andMe" ~ "Parkinson's disease\n(PD)",
phenotype == "BIP2021" ~ "Bipolar disorder (BIP)",
phenotype == "MDD2019" ~ "Major depressive disorder\n(MDD)",
phenotype == "SCZ2018" ~ "Schizophrenia\n(SCZ)"
)
)
plot <-
data_to_plot %>%
dplyr::mutate(
n_text_col =
case_when(
total > 200009 ~ "high",
TRUE ~ "low"
),
group = stringr::str_to_title(group)
) %>%
ggplot(
aes(
x = group,
y =
fct_relevel(
phenotype,
c(data_to_plot %>% dplyr::arrange(phenotype_group, phenotype) %>% dplyr::pull(phenotype) %>% unique() %>% rev())
),
fill = total
)
) +
geom_tile(colour = "black") +
geom_text(
aes(
label = total %>% scales::comma(),
colour = n_text_col
),
size = 3
) +
coord_equal() +
scale_x_discrete(
position = "top"
) +
labs(x = "", y = "") +
scale_fill_distiller(palette = "Greys", direction = 1) +
scale_colour_manual(values = c("white", "black")) +
guides(colour = F, fill = F) +
theme_minimal() +
theme(
panel.grid = element_blank(),
)
plot
Figure 4.1: Heatmap of samples numbers in each GWAS.
plot_global_corr(
global_corr = global,
n_phenotypes = 6
)
Figure 4.2: Global genetic correlations between pairs of phenotypes. Significant negative and positive correlations are indicated by blue and red fill, respectively. Non-significant correlations (p >= 0.05/n_tests) have a grey fill.
bivar_bonf %>%
dplyr::count(phen1) %>%
dplyr::bind_rows(
bivar_bonf %>%
dplyr::count(phen2) %>%
dplyr::rename(phen1 = phen2)
) %>%
dplyr::group_by(phen1) %>%
dplyr::summarise(
n_bivar = sum(n)
) %>%
dplyr::inner_join(
read_delim(
file = file.path(here::here("results", "01_input_prep"),
"input.info.txt"),
delim = "\t"
) %>%
dplyr::mutate(
phenotype = phenotype %>%
str_replace_all("[:digit:]", "") %>%
str_remove("\\..*"),
n_individ = cases + controls
),
by = c("phen1" = "phenotype")
) %>%
ggplot(
aes(
x = n_bivar,
y = n_individ,
colour = fct_relevel(
phen1,
fct_disease
)
)
) +
geom_point() +
scale_colour_brewer(
palette = "BrBG",
type = "div",
name = "GWAS"
) +
labs(
x = "Number of significant local genetic correlations",
y = "Sample size"
) +
theme_rhr
Figure 4.3: Plot of total number of significant bivariate local genetic correlations for each phenotype across all LD blocks versus GWAS sample size.
plot_bivar_chord_diagram(
bivar_corr = bivar_bonf,
fct_phen = fct_disease,
palette = RColorBrewer::brewer.pal(n = 6, name = "BrBG")
)
Figure 4.4: Chord diagram showing the number of significant bivariate local genetic correlations (p < 0.05/n_tests) between each of the phenotypes across all LD blocks. Positive and negative correlations are coloured red and blue, respectively.
4.3 Signficant phenotype correlations at LD block level
4.3.1 Text
- The LD block with the highest number of significant bivariate local \(r_{g}\)'s, LD block 758 and 1719, were positioned on chromosome 4 and 11, respectively (Figure 4.5, 4.6).
- The locus on chromosome 4 contains three protein-coding genes (GPM6A, WDR17 and SPATA4), while the LD block on chromosome 11 contains ten protein-coding genes, one of which includes DRD2, which encodes dopamine receptor D2. As perhaps expected, all three neuropsychiatric disorders positively correlate with one another at LD block 1719. Somewhat unexpected is the significant negative correlation between LBD and SCZ at LD block 1719.
print("Top LD block when Bonferroni-correcting (p < 0.05/n_bivar_tests):")## [1] "Top LD block when Bonferroni-correcting (p < 0.05/n_bivar_tests):"bivar_bonf %>%
dplyr::count(locus, chr, start, stop, n_snps) %>%
dplyr::arrange(-n) %>%
dplyr::inner_join(loci %>%
janitor::clean_names() %>%
dplyr::rename(
locus = loc
)) %>%
dplyr::slice_max(order_by = n, n = 1)4.3.2 Figures
results$bivar %>%
dplyr::filter(locus %in% unique(bivar_bonf$locus)) %>%
dplyr::mutate(
rho_fill =
case_when(
p < bivar_thres ~ round(rho, 2)
)
) %>%
ggplot(
aes(
x = phen1,
y = phen2,
fill = rho_fill,
label = round(rho, 2)
)
) +
geom_tile(colour = "black") +
geom_text(
size = 2
) +
facet_wrap(vars(locus), ncol = 5) +
scale_fill_distiller(palette = "RdBu", direction = -1, na.value = "#cccccc", limits = c(-1, 1)) +
theme_rhr
Figure 4.5: Heatmaps showing the standardised coefficient for genetic correlation (rho) for each association within the LD block with significant bivariate local genetic correlations (p < 0.05/n_tests). Significant negative and positive correlations are indicated by blue and red fill, respectively. Non-significant correlations have a grey fill.
plot_edge_diagram(
bivar_corr =
bivar_bonf,
phen = fct_disease,
ncol = 3,
)
Figure 4.6: Edge diagrams for each LD block showing the standardised coefficient for genetic correlation (rho) for each significant bivariate local genetic correlations (p < 0.05/n_tests). Significant negative and positive correlations are indicated by blue and red colour, respectively.
loci_gr <-
bivar_bonf %>%
dplyr::count(locus, chr, start, stop, n_snps) %>%
dplyr::arrange(locus) %>%
GenomicRanges::makeGRangesFromDataFrame(
.,
keep.extra.columns = TRUE,
ignore.strand = TRUE,
seqinfo = NULL,
seqnames.field = "chr",
start.field = "start",
end.field = "stop"
)
fig_list = vector(mode = "list", length = length(loci_gr))
for(i in 1:length(loci_gr)){
fig_list[[i]] <-
plot_locus(
locus_gr = loci_gr[i],
ref = ref
)
names(fig_list)[i] <- str_c("locus_", loci_gr[i]$locus)
}
fig_list## $locus_10
##
## $locus_63
##
## $locus_145
##
## $locus_266
##
## $locus_335
##
## $locus_375
##
## $locus_378
##
## $locus_385
##
## $locus_457
##
## $locus_472
##
## $locus_482
##
## $locus_483
##
## $locus_681
##
## $locus_711
##
## $locus_758
##
## $locus_891
##
## $locus_949
##
## $locus_950
##
## $locus_951
##
## $locus_952
##
## $locus_960
##
## $locus_1025
##
## $locus_1055
##
## $locus_1186
##
## $locus_1190
##
## $locus_1194
##
## $locus_1209
##
## $locus_1273
##
## $locus_1367
##
## $locus_1445
##
## $locus_1498
##
## $locus_1582
##
## $locus_1719
##
## $locus_1726
##
## $locus_1729
##
## $locus_1735
##
## $locus_1738
##
## $locus_1742
##
## $locus_1823
##
## $locus_1833
##
## $locus_1840
##
## $locus_1851
##
## $locus_1919
##
## $locus_1963
##
## $locus_2001
##
## $locus_2026
##
## $locus_2048
##
## $locus_2054
##
## $locus_2066
##
## $locus_2083
##
## $locus_2113
##
## $locus_2117
##
## $locus_2128
##
## $locus_2192
##
## $locus_2281
##
## $locus_2351
##
## $locus_2402
##
## $locus_2408
##
## $locus_2483
4.4 LD blocks with known neurodegenerative pleiotropic genes
4.4.1 Text
- Pleiotropic genes include:
- APOE (implicated in AD and LBD)
- TMEM175 (implicated in PD, LBD, and RBD)
- SNCA (implicated in PD, LBD and RBD)
- BIN1 (implicated in AD and LBD)
- GRN (implicated in AD, PD, FTD, and LBD)
- KIF5A (implicated in PD and ALS/FTD)
- AD and PD had sufficient local genetic signal (i.e. univariate heritability) at all tested LD blocks. LBD had insufficient local heritability in the GRN-, KIF5A and TMEM175-containing LD blocks, and therefore was not carried forward in these loci for bivariate testing.
- Of the LD blocks that contained implicated genes, only the SNCA- and APOE-containing LD blocks contained local \(r_{g}\)s that passed the Bonferroni significance cut-off (p < 0.05/n_tests; Figure 4.7a). Using a more lenient FDR correction, local \(r_{g}\)s were observed across all LD blocks, except the GRN- and TMEM175-containing LD blocks (Figure 4.7b).
- BIN1-containing LD block 319
- Bonferroni: no local \(r_{g}\)s.
- FDR: positive local \(r_{g}\) between AD and LBD.
- SNCA-containing LD block 681
- Bonferroni & FDR: negative local \(r_{g}\) between AD and PD.
- KIF5A-containing LD block 1794
- Bonferroni: no local \(r_{g}\)s.
- FDR: positive local \(r_{g}\)s between BIP and SCZ.
- Bonferroni: no local \(r_{g}\)s.
- APOE-containing LD block 2351
- Bonferroni: positive local \(r_{g}\) observed between AD and LBD and negative local \(r_{g}\) between PD and LBD.
- FDR: as above and a negative local \(r_{g}\) betyween AD and PD.
- BIN1-containing LD block 319
- In terms of expectations, APOE and BIN1 are what we would expect i.e. positive local \(r_{g}\) between AD and LBD, but some unexpected local \(r_{g}\)s too e.g. negative local \(r_{g}\) between AD and PD in APOE- and SNCA-containing LD blocks.
4.4.2 Tables
print("LD blocks containing pleiotropic genes:")## [1] "LD blocks containing pleiotropic genes:"loci_genes_filtered %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Univariate results filtered for neurodegenerative phenotypes AND LD blocks containing pleiotropic genes. Significant when p < 0.05/n_ld_blocks.")## [1] "Univariate results filtered for neurodegenerative phenotypes AND LD blocks containing pleiotropic genes. Significant when p < 0.05/n_ld_blocks."results$univ %>%
dplyr::inner_join(
loci_genes_filtered %>%
dplyr::select(locus, gene_name),
by = c("locus")
) %>%
dplyr::filter(
phen %in% c("AD", "LBD", "PD")
) %>%
dplyr::mutate(
significant =
case_when(
p < 0.05/nrow(loci) ~ TRUE,
TRUE ~ FALSE
)
) %>%
dplyr::select(
locus, gene_name, significant, everything()
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Bivariate results filtered for neurodegenerative phenotypes AND LD blocks containing pleiotropic genes.")## [1] "Bivariate results filtered for neurodegenerative phenotypes AND LD blocks containing pleiotropic genes."results$bivar %>%
dplyr::inner_join(
loci_genes_filtered %>%
dplyr::select(locus, gene_name),
by = c("locus")
) %>%
dplyr::filter(
phen1 %in% c("AD", "LBD", "PD") & phen2 %in% c("AD", "LBD", "PD")
) %>%
dplyr::mutate(
bonferroni =
case_when(
p < bivar_thres ~ TRUE,
TRUE ~ FALSE
),
fdr =
case_when(
fdr < 0.05 ~ TRUE,
TRUE ~ FALSE
)
) %>%
dplyr::select(
locus, gene_name, bonferroni, fdr, everything()
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.4.3 Figures
a <-
plot_edge_diagram(
bivar_corr =
results$bivar %>%
dplyr::inner_join(
loci_genes_filtered %>%
dplyr::select(locus, gene_name),
by = c("locus")
) %>%
dplyr::mutate(
locus = str_c(locus, ": contains ", gene_name)
),
p_threshold = bivar_thres,
multiple_corr = FALSE,
phen = fct_disease,
ncol = 3
)
b <-
plot_edge_diagram(
bivar_corr =
results$bivar %>%
dplyr::inner_join(
loci_genes_filtered %>%
dplyr::select(locus, gene_name),
by = c("locus")
) %>%
dplyr::mutate(
locus = str_c(locus, ": contains ", gene_name)
) %>%
dplyr::filter(fdr < 0.05),
multiple_corr = FALSE,
phen = fct_disease,
ncol = 3
) +
theme(legend.position = "none")
cowplot::plot_grid(
a,b,
ncol = 1,
rel_heights = c(1,2),
align = "v",
axis = "lr",
labels = letters[1:2]
)
Figure 4.7: Edge diagrams for each LD block showing local genetic correlations when corrected using (a) the stringent Bonferroni correction or (b) the lenient FDR correction. The standardised coefficient for genetic correlation (rho) is indicated for significant local correlations. Negative and positive correlations are indicated by blue and red colour, respectively.
5 Session info
## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.5 (2021-03-31)
## os Ubuntu 16.04.7 LTS
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate en_GB.UTF-8
## ctype en_GB.UTF-8
## tz Europe/London
## date 2022-02-09
##
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
## package * version date lib source
## assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.0.5)
## backports 1.3.0 2021-10-27 [1] CRAN (R 4.0.5)
## Biobase 2.50.0 2020-10-27 [2] Bioconductor
## BiocGenerics * 0.36.1 2021-04-16 [2] Bioconductor
## BiocParallel 1.24.1 2020-11-06 [2] Bioconductor
## Biostrings 2.58.0 2020-10-27 [2] Bioconductor
## bit 4.0.4 2020-08-04 [2] CRAN (R 4.0.5)
## bit64 4.0.5 2020-08-30 [2] CRAN (R 4.0.5)
## bitops 1.0-7 2021-04-24 [2] CRAN (R 4.0.5)
## bookdown 0.22 2021-04-22 [2] CRAN (R 4.0.5)
## broom 0.7.10 2021-10-31 [1] CRAN (R 4.0.5)
## bslib 0.3.1 2021-10-06 [1] CRAN (R 4.0.5)
## cellranger 1.1.0 2016-07-27 [2] CRAN (R 4.0.5)
## chron 2.3-56 2020-08-18 [1] CRAN (R 4.0.5)
## circlize * 0.4.13 2021-06-09 [1] CRAN (R 4.0.5)
## cli 3.1.0 2021-10-27 [1] CRAN (R 4.0.5)
## colorspace 2.0-2 2021-06-24 [2] CRAN (R 4.0.5)
## cowplot 1.1.1 2020-12-30 [2] CRAN (R 4.0.5)
## crayon 1.4.2 2021-10-29 [1] CRAN (R 4.0.5)
## crosstalk 1.1.1 2021-01-12 [2] CRAN (R 4.0.5)
## data.table 1.14.2 2021-09-27 [1] CRAN (R 4.0.5)
## DBI 1.1.1 2021-01-15 [2] CRAN (R 4.0.5)
## dbplyr 2.1.1 2021-04-06 [2] CRAN (R 4.0.5)
## DelayedArray 0.16.3 2021-03-24 [2] Bioconductor
## digest 0.6.29 2021-12-01 [1] CRAN (R 4.0.5)
## dplyr * 1.0.7 2021-06-18 [2] CRAN (R 4.0.5)
## DT 0.19 2021-09-02 [1] CRAN (R 4.0.5)
## ellipsis 0.3.2 2021-04-29 [2] CRAN (R 4.0.5)
## evaluate 0.14 2019-05-28 [2] CRAN (R 4.0.5)
## fansi 0.5.0 2021-05-25 [2] CRAN (R 4.0.5)
## farver 2.1.0 2021-02-28 [2] CRAN (R 4.0.5)
## fastmap 1.1.0 2021-01-25 [2] CRAN (R 4.0.5)
## forcats * 0.5.1 2021-01-27 [2] CRAN (R 4.0.5)
## fs 1.5.1 2021-11-30 [1] CRAN (R 4.0.5)
## generics 0.1.1 2021-10-25 [1] CRAN (R 4.0.5)
## GenomeInfoDb * 1.26.7 2021-04-08 [2] Bioconductor
## GenomeInfoDbData 1.2.4 2021-07-03 [2] Bioconductor
## GenomicAlignments 1.26.0 2020-10-27 [2] Bioconductor
## GenomicRanges * 1.42.0 2020-10-27 [2] Bioconductor
## ggforce 0.3.3 2021-03-05 [2] CRAN (R 4.0.5)
## ggplot2 * 3.3.5 2021-06-25 [2] CRAN (R 4.0.5)
## ggraph * 2.0.5 2021-02-23 [1] CRAN (R 4.0.5)
## ggrepel 0.9.1 2021-01-15 [2] CRAN (R 4.0.5)
## GlobalOptions 0.1.2 2020-06-10 [1] CRAN (R 4.0.5)
## glue 1.5.1 2021-11-30 [1] CRAN (R 4.0.5)
## graphlayouts 0.7.1 2020-10-26 [1] CRAN (R 4.0.5)
## gridExtra 2.3 2017-09-09 [2] CRAN (R 4.0.5)
## gtable 0.3.0 2019-03-25 [2] CRAN (R 4.0.5)
## haven 2.4.3 2021-08-04 [1] CRAN (R 4.0.5)
## here 1.0.1 2020-12-13 [1] CRAN (R 4.0.5)
## highr 0.9 2021-04-16 [2] CRAN (R 4.0.5)
## hms 1.1.1 2021-09-26 [1] CRAN (R 4.0.5)
## htmltools 0.5.2 2021-08-25 [1] CRAN (R 4.0.5)
## htmlwidgets 1.5.4 2021-09-08 [1] CRAN (R 4.0.5)
## httr 1.4.2 2020-07-20 [2] CRAN (R 4.0.5)
## igraph 1.2.7 2021-10-15 [1] CRAN (R 4.0.5)
## IRanges * 2.24.1 2020-12-12 [2] Bioconductor
## janitor * 2.1.0 2021-01-05 [1] CRAN (R 4.0.5)
## jquerylib 0.1.4 2021-04-26 [2] CRAN (R 4.0.5)
## jsonlite 1.7.2 2020-12-09 [2] CRAN (R 4.0.5)
## knitr 1.36 2021-09-29 [1] CRAN (R 4.0.5)
## labeling 0.4.2 2020-10-20 [2] CRAN (R 4.0.5)
## lattice 0.20-44 2021-05-02 [2] CRAN (R 4.0.5)
## LAVA 0.0.6 2021-09-01 [1] xgit (@7be3421)
## lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.0.5)
## lubridate 1.8.0 2021-10-07 [1] CRAN (R 4.0.5)
## magrittr 2.0.1 2020-11-17 [2] CRAN (R 4.0.5)
## MASS 7.3-54 2021-05-03 [2] CRAN (R 4.0.5)
## Matrix 1.3-4 2021-06-01 [2] CRAN (R 4.0.5)
## MatrixGenerics 1.2.1 2021-01-30 [2] Bioconductor
## matrixStats 0.61.0 2021-09-17 [1] CRAN (R 4.0.5)
## modelr 0.1.8 2020-05-19 [2] CRAN (R 4.0.5)
## munsell 0.5.0 2018-06-12 [2] CRAN (R 4.0.5)
## pillar 1.6.4 2021-10-18 [1] CRAN (R 4.0.5)
## pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.0.5)
## polyclip 1.10-0 2019-03-14 [2] CRAN (R 4.0.5)
## purrr * 0.3.4 2020-04-17 [2] CRAN (R 4.0.5)
## qdapTools 1.3.5 2020-04-17 [1] CRAN (R 4.0.5)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.0.5)
## RColorBrewer 1.1-2 2014-12-07 [2] CRAN (R 4.0.5)
## Rcpp 1.0.7 2021-07-07 [1] CRAN (R 4.0.5)
## RCurl 1.98-1.5 2021-09-17 [1] CRAN (R 4.0.5)
## readr * 2.0.2 2021-09-27 [1] CRAN (R 4.0.5)
## readxl 1.3.1 2019-03-13 [2] CRAN (R 4.0.5)
## reprex 2.0.0 2021-04-02 [2] CRAN (R 4.0.5)
## rlang 0.4.12 2021-10-18 [1] CRAN (R 4.0.5)
## rmarkdown 2.11 2021-09-14 [1] CRAN (R 4.0.5)
## rprojroot 2.0.2 2020-11-15 [2] CRAN (R 4.0.5)
## Rsamtools 2.6.0 2020-10-27 [2] Bioconductor
## rstudioapi 0.13 2020-11-12 [2] CRAN (R 4.0.5)
## rtracklayer * 1.50.0 2020-10-27 [2] Bioconductor
## rvest 1.0.1 2021-07-26 [2] CRAN (R 4.0.5)
## S4Vectors * 0.28.1 2020-12-09 [2] Bioconductor
## sass 0.4.0 2021-05-12 [2] CRAN (R 4.0.5)
## scales 1.1.1 2020-05-11 [2] CRAN (R 4.0.5)
## sessioninfo * 1.1.1 2018-11-05 [2] CRAN (R 4.0.5)
## shape 1.4.6 2021-05-19 [1] CRAN (R 4.0.5)
## snakecase 0.11.0 2019-05-25 [1] CRAN (R 4.0.5)
## stringi 1.7.6 2021-11-29 [1] CRAN (R 4.0.5)
## stringr * 1.4.0 2019-02-10 [2] CRAN (R 4.0.5)
## SummarizedExperiment 1.20.0 2020-10-27 [2] Bioconductor
## tibble * 3.1.6 2021-11-07 [1] CRAN (R 4.0.5)
## tidygraph 1.2.0 2020-05-12 [1] CRAN (R 4.0.5)
## tidyr * 1.1.4 2021-09-27 [1] CRAN (R 4.0.5)
## tidyselect 1.1.1 2021-04-30 [2] CRAN (R 4.0.5)
## tidyverse * 1.3.1 2021-04-15 [1] CRAN (R 4.0.5)
## tweenr 1.0.2 2021-03-23 [2] CRAN (R 4.0.5)
## tzdb 0.2.0 2021-10-27 [1] CRAN (R 4.0.5)
## utf8 1.2.2 2021-07-24 [2] CRAN (R 4.0.5)
## vctrs 0.3.8 2021-04-29 [2] CRAN (R 4.0.5)
## viridis 0.6.2 2021-10-13 [1] CRAN (R 4.0.5)
## viridisLite 0.4.0 2021-04-13 [2] CRAN (R 4.0.5)
## vroom 1.5.5 2021-09-14 [1] CRAN (R 4.0.5)
## withr 2.4.3 2021-11-30 [1] CRAN (R 4.0.5)
## xfun 0.27 2021-10-18 [1] CRAN (R 4.0.5)
## XML 3.99-0.8 2021-09-17 [1] CRAN (R 4.0.5)
## xml2 1.3.2 2020-04-23 [2] CRAN (R 4.0.5)
## XVector 0.30.0 2020-10-27 [2] Bioconductor
## yaml 2.2.1 2020-02-01 [2] CRAN (R 4.0.5)
## zlibbioc 1.36.0 2020-10-27 [2] Bioconductor
##
## [1] /home/rreynolds/R/x86_64-pc-linux-gnu-library/4.0
## [2] /opt/R/4.0.5/lib/R/library